import numpy as np  
import matplotlib.pyplot as plt  
from scipy.stats import norm, gaussian_kde  
  
# 1. 生成x和y的随机样本  
np.random.seed(0)  # 为了结果的可重复性  
n_samples = 1000000  # 样本数量  
x = np.random.normal(0, 1, n_samples)  
y = np.random.normal(0, 1, n_samples)  
  
# 排除y为0的情况（因为z = x / y在y=0时无定义）  
x_filtered = x[y != 0]  
y_filtered = y[y != 0]  
z = x_filtered * y_filtered  
  
# 2. 使用核密度估计来近似z的PDF  
z_kde = gaussian_kde(z)  
x_values_z = np.linspace(min(z), max(z), 10000)  
z_pdf = z_kde.evaluate(x_values_z)  
  
# 3. 使用直方图来估计z的PDF  
bins = np.linspace(min(z), max(z), 30000)  # 选择合适的bin范围和数量  
plt.hist(z, bins=bins, density=True, alpha=0.6, label='z = x / y')  
  
# 4. 绘制x和y的PDF图像（作为参考）  
x_values = np.linspace(-4, 4, 10000)  
plt.plot(x_values, norm.pdf(x_values), 'r', lw=2, label='x (standard normal)')  
plt.plot(x_values, norm.pdf(x_values), 'g', lw=2, label='y (standard normal)')  # 这里为了清晰起见，通常不会同时绘制  
  
# 5. 绘制z的PDF近似  
plt.plot(x_values_z, z_pdf, 'b', lw=2, label='z = x / y (KDE)')  
  
# 6. 设置图例、标题等  
plt.legend()  
plt.title('PDFs of x, y, and z = x / y')  
plt.xlabel('Value')  
plt.ylabel('Probability Density')  
plt.show()